# Replication files for:
# "Generalizability of Heterogeneous Treatment Effect Estimates Across Samples"
# Alexander Coppock, Thomas J. Leeper, and Kevin J. Mullinix
# Proceedings of the National Academy of Sciences, Forthcoming

rm(list = ls())
library(dplyr)
library(tidyr)
library(readr)
library(estimatr)
library(lmtest)

# Superordinate IDs -------------------------------------------------------

superordinate_id_stacked <- read_rds("superordinate_id_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (Z_particularism + age_3 + pid_3) + sample * (Z_particularism + age_3 + pid_3), data = superordinate_id_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (Z_particularism + age_3 + pid_3) * sample, data = superordinate_id_stacked)

superordinate_id <- waldtest(fit_1, fit_2, test = "F")

# Patriot Act -------------------------------------------------------------

patriot_act_stacked <- read_rds("patriot_act_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female), data = patriot_act_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female) * sample, data = patriot_act_stacked)

patriot_act <- waldtest(fit_1, fit_2, test = "F")

# Elite Endorsements ------------------------------------------------------

elite_endorsements_stacked <- read_rds("elite_endorsements_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female), data = filter(elite_endorsements_stacked, Z_imm_match != "control"))
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female)* sample, data = filter(elite_endorsements_stacked, Z_imm_match != "control"))

elite_endorsements <- waldtest(fit_1, fit_2, test = "F")

# Free Trade --------------------------------------------------------------

free_trade_stacked <- read_rds("free_trade_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) +
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = free_trade_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = free_trade_stacked)

free_trade <- waldtest(fit_1, fit_2, test = "F")

# Frame Breadth -----------------------------------------------------------

frame_breadth_stacked <- read_rds("frame_breadth_stacked.rds")
frame_breadth_stacked <- frame_breadth_stacked %>% filter(sample != "gfk")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = frame_breadth_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = frame_breadth_stacked)

frame_breadth <- waldtest(fit_1, fit_2, test = "F")

# Polarization ------------------------------------------------------------

polarization_stacked <- read_rds("polarization_stacked.rds")

polarization_stacked <- polarization_stacked %>% filter(sample != "gfk")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) +
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = polarization_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = polarization_stacked)

polarization <- waldtest(fit_1, fit_2, test = "F")

# Immigration -------------------------------------------------------------

immigration_stacked <- read_rds("immigration_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + educ_3 + female) + 
                     sample * (age_3 + pid_3 + educ_3 + female), data = immigration_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + educ_3 + female) * sample, data = immigration_stacked)

immigration <- waldtest(fit_1, fit_2, test = "F")

# System Threat -----------------------------------------------------------

system_threat_stacked <- read_rds("system_threat_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = system_threat_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = system_threat_stacked)

system_threat <- waldtest(fit_1, fit_2, test = "F")

# Expert Economists -------------------------------------------------------

expert_economists_stacked_long <- read_rds("expert_economists_stacked_long.rds")

expert_economists_stacked_long <- 
  expert_economists_stacked_long %>%
  filter(!is.na(Z), !is.na(Y_s))
  
fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), clusters = anon_ID, data = expert_economists_stacked_long)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, clusters = anon_ID, data = expert_economists_stacked_long)

expert_economists <- waldtest(fit_1, fit_2, test = "F")
  
# Mental Illness ----------------------------------------------------------

mental_illness_stacked <- read_rds("mental_illness_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = mental_illness_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = mental_illness_stacked)

mental_illness <- waldtest(fit_1, fit_2, test = "F")

# Death Penalty -----------------------------------------------------------

death_penalty_stacked <- read_rds("death_penalty_stacked.rds")
death_penalty_stacked <- 
  death_penalty_stacked %>%
  filter(Z_CP_3 != "control")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = death_penalty_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = death_penalty_stacked)

death_penalty <- waldtest(fit_1, fit_2, test = "F")

# berganS20 ---------------------------------------------------------------

berganS20_stacked <- read_rds("berganS20_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = berganS20_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = berganS20_stacked)

berganS20 <- waldtest(fit_1, fit_2, test = "F")


# brandtS1 ----------------------------------------------------------------

brandtS1_stacked <- read_rds("brandtS1_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female), data = brandtS1_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female) * sample, data = brandtS1_stacked)

brandtS1 <- waldtest(fit_1, fit_2, test = "F")


# caprarielloS2 -----------------------------------------------------------

caprarielloS2_stacked <- read_rds("caprarielloS2_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = caprarielloS2_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = caprarielloS2_stacked)

caprarielloS2 <- waldtest(fit_1, fit_2, test = "F")

# converseS16 -------------------------------------------------------------

converseS16_stacked <- read_rds("converseS16_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + white + educ_3 + female) + 
                     sample * (age_3 + white + educ_3 + female), data = converseS16_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + white + educ_3 + female) * sample, data = converseS16_stacked)

converseS16 <- waldtest(fit_1, fit_2, test = "F")

# dennyS17 ----------------------------------------------------------------

dennyS17_stacked <- read_rds("dennyS17_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = dennyS17_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = dennyS17_stacked)

dennyS17 <- waldtest(fit_1, fit_2, test = "F")

# flavinS4 ----------------------------------------------------------------

flavinS4_stacked <- read_rds("flavinS4_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = flavinS4_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = flavinS4_stacked)

flavinS4 <- waldtest(fit_1, fit_2, test = "F")

# gashS5 ------------------------------------------------------------------

gashS5_stacked <- read_rds("gashS5_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = gashS5_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = gashS5_stacked)

gashS5 <- waldtest(fit_1, fit_2, test = "F")


# jacobsenS7 --------------------------------------------------------------

jacobsenS7_stacked <- read_rds("jacobsenS7_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = jacobsenS7_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = jacobsenS7_stacked)

jacobsenS7 <- waldtest(fit_1, fit_2, test = "F")


# melloS6 -----------------------------------------------------------------

melloS6_stacked <- read_rds("melloS6_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + white + educ_3 + female) + 
                     sample * (age_3 + white + educ_3 + female), data = melloS6_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + white + educ_3 + female) * sample, data = melloS6_stacked)

melloS6 <- waldtest(fit_1, fit_2, test = "F")

# parmerS15 ---------------------------------------------------------------

parmerS15_stacked <- read_rds("parmerS15_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = parmerS15_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = parmerS15_stacked)

parmerS15 <- waldtest(fit_1, fit_2, test = "F")

# pedullaS18 --------------------------------------------------------------
pedullaS18_stacked <- read_rds("pedullaS18_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = pedullaS18_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = pedullaS18_stacked)

pedullaS18 <- waldtest(fit_1, fit_2, test = "F")

# piazzaS8 ----------------------------------------------------------------
piazzaS8_stacked <- read_rds("piazzaS8_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = piazzaS8_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = piazzaS8_stacked)

piazzaS8 <- waldtest(fit_1, fit_2, test = "F")

# shaferS9 ----------------------------------------------------------------
shaferS9_stacked <- read_rds("shaferS9_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = shaferS9_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = shaferS9_stacked)

shaferS9 <- waldtest(fit_1, fit_2, test = "F")

# thompsonS10 -------------------------------------------------------------
thompsonS10_stacked <- read_rds("thompsonS10_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = thompsonS10_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = thompsonS10_stacked)

thompsonS10 <- waldtest(fit_1, fit_2, test = "F")

# turagaS11 ---------------------------------------------------------------
turagaS11_stacked <- read_rds("turagaS11_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = turagaS11_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = turagaS11_stacked)

turagaS11 <- waldtest(fit_1, fit_2, test = "F")

# wallaceS12 --------------------------------------------------------------

wallaceS12_stacked <- read_rds("wallaceS12_stacked.rds")

fit_1 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) + 
                     sample * (age_3 + pid_3 + white + educ_3 + female + ideo_3), data = wallaceS12_stacked)
fit_2 <- lm_robust(Y_s ~ Z * (age_3 + pid_3 + white + educ_3 + female + ideo_3) * sample, data = wallaceS12_stacked)

wallaceS12 <- waldtest(fit_1, fit_2, test = "F")


# bind and merge ----------------------------------------------------------

get_f <- function(x){data_frame(p_value = x$`Pr(>F)`[2])}

f_tests_df <-
  bind_rows(
    superordinate_id = get_f(superordinate_id),
    patriot_act = get_f(patriot_act),
    elite_endorsements = get_f(elite_endorsements),
    free_trade = get_f(free_trade),
    frame_breadth = get_f(frame_breadth),
    polarization = get_f(polarization),
    immigration = get_f(immigration),
    system_threat = get_f(system_threat),
    expert_economists = get_f(expert_economists),
    mental_illness = get_f(mental_illness),
    death_penalty = get_f(death_penalty),
    berganS20 = get_f(berganS20),
    brandtS1 = get_f(brandtS1),
    caprarielloS2 = get_f(caprarielloS2),
    converseS16 = get_f(converseS16),
    dennyS17 = get_f(dennyS17),
    flavinS4 = get_f(flavinS4),
    gashS5 = get_f(gashS5),
    jacobsenS7 = get_f(jacobsenS7),
    melloS6 = get_f(melloS6),
    parmerS15 = get_f(parmerS15),
    pedullaS18 = get_f(pedullaS18),
    piazzaS8 = get_f(piazzaS8),
    shaferS9 = get_f(shaferS9),
    thompsonS10 = get_f(thompsonS10),
    turagaS11 = get_f(turagaS11),
    wallaceS12 = get_f(wallaceS12),
    .id = "study"
  )


write_rds(f_tests_df, path = "CLM_f_tests_df.rds")

# Consistent with Figure 2, we fail to reject the null hypothesis most of the time 
# (25 of 27 opportunities), indicating that whatever heterogeneity in treatment effects there may be,
# the patterns do not differ greatly across samples.
f_tests_df %>% summarize(sum(p_value > 0.05), n = n())




